Journal of the Bombay Natural History Society, 106(1), Jan-Apr 2009 


38-44 


GENETIC DIFFERENTIATION OF ARGALI SHEEP OVIS AMMON IN MONGOLIA 
REVEALED BY MITOCHONDRIAL CONTROL REGION 
AND NUCLEAR MICROSATELLITES ANALYSES 


Jru Fenc!, MICHAEL R. FRISINA?, MICHAEL S. WEBSTER? AND GOMBOSUREN ULZzIIMAA* 


'Department of Biological Sciences, 109 Cooke Hall, State University of New York at Buffalo, Buffalo, NY 14260, USA. 
*Montana Department of Fish, Wildlife and Parks, 1330 West Gold Street, Butte, MT 59701-2112, USA. 


Email: habitat@bresnan.net 


‘Department of Biological Sciences, 109 Cooke Hall, State University of New York at Buffalo, NY 14260, 


Washington State University, Pullman, WA 99164-4236, USA. 


‘Mongolian National Agriculture University, Ulaanbaatar, Mongolia. 


The genetic distinctiveness and possible gene flow among the Argali Sheep (Ovis ammon) populations in Mongolia have 
been controversial, due to a high degree of morphological variation among populations and an apparent lack of physical 
barriers to dispersal. We studied the population genetic structure of Argali sheep in Mongolia using both mitochondrial 
control region sequences (613 bp) and 14 nuclear microsatellite markers. Mitochondrial results suggest two evolutionarily 
distinct lineages, one in the Altay Moun‘ains and the other in the Hangay Mountains and eastern Gobi Desert. 
Microsatellite analysis indicated genetic differentiation among these three regions, and also indicated similar levels of 
genetic differentiation and gene flow among all pair-wise comparisons. These results suggest genetic differentiation 
among the Mongolian populations of this endangered mammal. 
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INTRODUCTION 


The Argali Sheep (Ovis ammon), which is the largest 
species of wild sheep in the world, has become endangered 
due to poaching and habitat destruction (Valdez 1982; U.S. 
Fish and Wildlife Service 1996). Currently in Mongolia, Argali 
Sheep are patchily distributed in the Altay Mountains, Hangay 
Mountains, and Gobi Desert (Mallon 1985). Following a 
country-wide survey in 2002, Frisina et al. (2007) estimated a 
Mongolian Argali population of about 20,000. Across this 
range, both elevation and habitat productivity decrease 
gradually without obvious physical barriers to dispersal 
(Frisina 1998). Morphologically, Mongolian Argali Sheep is 
highly variable, and there is a general trend for average body 
size to decrease as elevation decreases from west to east, 
with Altay Argali being the largest in body size of all O. ammon 
(Geist 1991). 

The variable morphology and the lack of obvious 
geographic barriers to gene flow have led to a long-debated 
controversy regarding the taxonomic status of Mongolian 
Argali and the delineation of genetically distinct populations. 
Allen (1940) considered all Mongolian Argali to be one 
subspecies O.a. darwini, but currently two subspecies are 
commonly recognized: O.a. ammon (Altay Argali) are large 
argalis from the Altay mountain region, and O.a. darwini (Gobi 
Argali) are smaller argalis from Gobi desert region (Sopin 1982; 
Valdez 1982; Geist 1991; Mitchell and Frisina 2007). Detailed 
analysis of cranial morphology show that O.a. ammon and 
O.a. darwini are morphologically distinct from other argalis 


and from each other, thus supporting subspecific recognition 
of these taxa (Kapitanova ef al, 2004). The taxonomic position 
of argali from the Hangay region of Mongolia is unclear: Sopin 
(1982) and Geist (1991) considered these argalis to be similar 
to those from Altay (O.a. ammon), but genetic analyses by 
Tserenbataa et al. (2004) suggest that they may be more closely 
allied to Gobi argalts (see below). 

Although a species of conservation concern, at present 
little is known about population structure and gene flow 
among argali populations in Mongolia. Tserenbataa et al. 
(2004) examined genetic variation at the mitochondrial ND5 
locus, and found little genetic differentiation among 
populations in Mongolia and nearby regions of Kazakhstan 
and Kyrgyzstan. What little genetic differentiation they did 
find appeared to separate Gobi/Hangay from Altay/ 
Kazakhstan/Kyrgyszstan. Tserenbataa ef al. (2004) attributed 
the lack of genetic differentiation to high levels of female- 
mediated gene flow among populations, and concluded that 
argali populations from all of Mongolia and nearby regions 
of China and Russia should be considered a single 
“evolutionary significant unit” (or subspecies) with two 
management units. 

The conclusions of Tserenbataa et al. (2004) contrasts 
with those from morphometric analyses, which suggest two 
argali subspecies in Mongolia (Kapitanova ef al. 2004). 
Moreover, the shallow ND5 phylogenetic tree presented by 
Tserenbataa et al. (2004) suggests that the lack of genetic 
differentiation among populations may be due to incomplete 
lineage sorting rather than to the movement of individuals 
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Fig. 1: Geographic distribution of sampling locations (filled squares) in Mongolia 


(Avise 2000). To help resolve this taxonomic controversy and 
aid in conservation efforts, we studied two molecular genetic 
markers that show relatively high evolutionary rates — the 
mitochondrial control region and nuclear microsatellites — 
for argali samples collected in Mongolia. 


Abbreviations: mtDNA = mitochondrial DNA, 
MP = maximum parsimony, ML = maximum likelihood, 


ME = minimum evolution 
MATERIALAND METHODS 


We collected a total of 58 argali samples from Altay 
Mountains, Hangay Mountains and eastern Gobi Desert in 
Mongolia (Fig. 1), including tissue samples (skin and liver) 
collected from legally hunted individuals and bone samples 
(horn fragments and teeth) collected from carcasses found in 
the field. Skin samples of Snow Sheep (Ovis nivicola) from 
Russia were collected to serve as an outgroup for phylogenetic 
comparisons. Despite considerable and repeated efforts to 
extract and amplify DNA from all sources, some sources 
(e.g., some bone samples from carcasses) proved difficult 
and did not yield usable DNA. Consequently, the sample 
sizes for mtDNA and microsatellite analyses (see below) differ 
from each other and from the total number of samples collected. 

We used a standard proteinase K digestion and phenol/ 
chloroform methods to extract genomic DNA (Sambrook ef 
al. 1989). The mtDNA control region was amplified via 
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polymerase chain reaction (PCR) using primers modified for 
ungulates (Murray et al. 1995), and sequenced the portion 
proximal to tRNA®® via cycle sequencing (Feng et al. 2001). 
Replicate amplifications were sequenced for most samples, 
and replicates always yielded identical results. Moreover, all 
sequences were clean and easily scored, suggesting that we 
did not co-amplify nuclear paralogs or encounter 
heteroplasmy. All haplotype sequences have been deposited 
in Genbank (accession numbers AY315886-AY215899). We 
used maximum parsimony, maximum likelihood, and 
minimum evolution approaches for phylogenetic analyses of 
the control region haplotypes (details given below). These 
analyses were conducted using PAUP*4.0b2 (Swofford 
1998). 

We screened 37 pairs of primers of dinucleotide-repeat 
microsatellites developed from domestic sheep (Crawford ef 
al. 1995) in 10 argali individuals, and found 14 loci to be 
polymorphic: ILS5, ILS56, MAF33, MAF36, MAF48, 
MAF64, MAF209, FCB128, FCB226, FCB304, OHH35, 
OHHS6, OVH72, and OVH116. Published primers and 
annealing temperatures (Crawford ef al. 1995) to amplify 
alleles under the following conditions: 10 ul total reaction 
volume containing 20-50 ug genomic DNA, 1x PCR buffer 
(Roche), 0.05 mM dNTPs, 0.05 mM of each primer, 3.0 mM 
MgCl, 0.5U Taq polymerase, and 15 wCi *P-dATP (to label 
alleles). The PCR program was 94 °C for 3 min, followed by 
35 cycles of 94 °C for 1 min, annealing temperature for | min, 
and 72 °C for 45 sec. The last cycle was followed by a 5 min 
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extension at 72 °C. We separated PCR products by 
electrophoresis through a 6% polyacrylamide gel and ran a 
M13 sequence standard along with our PCR products for 
sizing the alleles. As with sequence analyses, replicates were 
amplified and electrophoresed for most samples to ensure 
accuracy of results. 

For each population, the program GENEPOP web 
version 3.lc (Raymond and Rousset 1995a) was used to 
calculate genetic diversity for each locus, as well as the mean 
observed heterozygosity across all loci. Global tests of both 
allelic and genotypic distributions were performed to detect 
population differentiation. F , (calculated using FSTAT 1.2, 
Goudet 1995) rather than R_ was used to measure population 
differentiation because the former performs better when 
sample size is small (Gaggiotti et al. 1999). The genetic 
distance calculator at http://www.biology.ualberta.ca was also 
used to calculate pair-wise Nei’s genetic distances between 
populations. At the individual level, we conducted assignment 
tests (Paetkau et al. 1995; Waser and Strobeck 1998) using 
the calculator available at http://www.biology.ualberta.ca. We 
also calculated the pair-wise allele-sharing genetic distance 
(Bowcock et al. 1994) matrix, which was then subjected to 
PAUP*4.0b2 (Swofford 1998) and multidimensional scaling 
analysis in two dimensions (Manly 1997) with SPSS to test 
whether genetic similarity reflects geographic groupings. 


RESULTS 


Mitochondrial Control Region Phylogeny: We 
obtained 14 argali control region haplotypes from 
17 sequences. We aligned 613 bp (including indels), of which 
92 (15.0%) were variable (Fig. 2), and 32 of these were parsimony 


informative. The maximum likelihood estimate of transition/ 
transversion ratio was 3.4:1. Base frequencies did not differ 
significantly across taxa (%? = 23.07, df = 45, p = 0.997), with 
A=38.9%, C=21.8%, G= 11.3%, and T = 28.0%. The sequence 
data contained significant phylogenetic signal as indicated 
by both a permutation test (PTP test, 1000 replicates, 
p < 0.001), and a tree length skewness test (g/ test, 
10,000 random trees, p < 0.01). Hierarchical likelihood ratio 
tests indicated that the optimal sequence evolution model for 
our observed data was the HK Y+G model, which incorporates 
unequal base frequencies, unequal transition vs. transversion 
rates, and among site rate heterogeneity. The rate 
heterogeneity distribution parameter was & = 0.56 (S.E. = 0.10), 
and the total heterogeneity (Gu etal. 1995) was 0.64. 

Four maximum parsimony (MP) trees were found through 
a branch and bound search. Maximum likelihood (ML) and 
minimum evolution (ME) analyses yielded tree topologies 
that were concordant with the strict consensus MP tree (Fig. 
3). The tree topology indicated that haplotypes from Hangay 
and east Gobi are more closely related to each other than they 
are to those from Altay. Seven out of eight haplotypes from 
Altay formed a single well-supported clade (“Altay group”). 
Haplotypes from Hangay and Gobi, plus a single haplotype 
from southern Altay, formed another well-supported clade 
(““Hangay/Gobi group’). Pair-wise ML genetic distance 
between the Altay group and Hangay/Gobi group (5.32% 
+).08%) was greater than that within the Altay (1.19% 40.13%) 
and Hangay/Gobi groups (0.67% +0.07%). Factoring out intra- 
group variation, the average ML genetic distance between 
the two groups was 4.39%. 

Microsatellite Diversity and Differentiation: We 
obtained microsatellite genotype data at 14 loci for a total of 
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7M(A)  ATTTCGCTGCTCACATAACAACCCATACAGAAAAGCACAATCACTTAGGATGTCAATCGTTAC-TAACCCAGTAAAGTATAG-CATTTACCC 
PAS CR) oo CC A a cs a oe ee E C.GTCCGAAGCACTG.CT-G.GTT..... Peau TAC. ee A..A. 
M6 (A) CO PEATE M is ofA sea uss ae C.GTCCGAAG.ACTG.C...-GTT...T....... TAC ceuey eee ae ne 
M7 (A) FG mene oa i oae e e A a eee C.GTCCGAAG.ACTG.CT-G.GTT..... T...-.TAC......... -..A. 
M9 (A) BEC yn AO a e EA A A a ee C.GTCCGAAG.ACTG.CT-G.GTT......... =, TAC.. AG... A.. A. 
M19(A) .CC.TA..-------------------------------- C.GTCCGAAG.ACTG.CT-G.G.T....TT..... TAC. Vera -T.A. 
MOOI ie CCE A oa) 2 5. tay ee ae ae TGTCCGAA..ACTG.CT-G.GTTGT...T..... TAC ne AGAA. 
2 dA eee a E a a A C.GICCGAA. .ACTG.CT-G.GTT..... Peon PAC CT eee A..-. 
IS CY ere ee re OG acer a eng ere €.C..GC.-G. .TT.=G.. 2h. TA. C CAA 
MIA e e oea A E a a e a a nei ane cera CC GCO TI G. CN TACT CATA 
MIS (Ba T Eora A EREE e a eee C.C.. GC.CG.: TT.. G... A.-. TA. C CRCC ATA 
KIT C hea AA T E A A A E A C.C..GC.-G..TT.-G....AG..TA..C..... CGG =. 
MIS (eee Cissus occa se i a E a CC. .CC. -C -TT G- A: TA TCh. CAT 
Aro EE te ace de OY Saree, i E N C C. CC. G TI. C ACCT o CIAA. 


Fig. 2: Sequence alignment of polymorphic nucleotide sites for Argali control region haplotypes 
Note: Numbers along the top refer to positions in the consensus sequence; Sample names are given in the extreme left column; and 
letters in parentheses refer to sampling location (A = Altay, H = Hangay, and G = Gobi); the top row gives the sequence for a reference 
sample (7M); dots indicate nucleotides that are identical to the reference; and dashes indicate deletions 
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Fig. 3: Strict consensus tree of four most parsimony trees (tree 
length 343, Cl = 0.904, RI =0.907, RC = 0.820) rooted with snow 
sheep (Ovis nivicola). Transversion: transition ratio was 3:1, 
gaps were treated as a fifth state, and the branch and bound 
search algorithm was used. Numbers above the branches indicate 
the bootstrap values obtained through 1000 replications (only 
values above 50% are shown). The maximum likelinood tree had an 
identical topology, and numbers below branches indicate quartet 
puzzling support values (only values above 50% are shown) 


30 individuals (12 from Altay, 11 from Hangay, and 7 from 
Gobi). There were no significant deviations from Hardy- 
Weinberg equilibrium (U tests, Raymond and Rousset 1995b), 
and no significant linkage disequilibrium (Fisher exact tests). 

The total number of alleles for a locus across 
populations ranged from 4 to 13 (mean = 7.5, S.E. = 0.79), 
and observed heterozygosity ranged from 0.30 to 
0.77 (mean = 0.61, S.E. = 0.04). For Altay, Hangay, and Gobi 
populations, the average number of alleles per locus (+ S.E.) 
was 5.9 (+0.47), 5.0 (+0.47), and 3.6 (+0.36) respectively, 
and the mean observed heterozygosity was 0.61 (+ 0.04), 
0.65 (+ 0.06), and 0.53 (+ 0.07) respectively. 

Fisher’s exact test conducted on microsatellite allele 
and genotype frequencies showed significant population 
differentiation (P<0.001). The mean F for all loci was 
0.056 (S.D. =0.017), which was significantly greater than zero 
(permutation test with 1000 replications, P < 0.001). Pair-wise 
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F , values were similar across population pairs (range: 0.051 
to 0.056), as were Nei’s pair-wise genetic distances 
(range: 0.213 to 0.272). Likelihood assignment tests yielded a 
high percentage of correct assignments (26 of 30, 87%). 
For Gobi, all of 7 individuals were assigned to Gobi. For Altay, 
10 of 12 individuals were assigned to Altay, | was assigned to 
Hangay, and 1 was assigned to Gobi. For Hangay, 9 of 
11 individuals were assigned to Hangay, | was assigned to 
Altay, and | was assigned to Gobi. The percentage of correct 
assignment was significantly higher than expected by chance 
(y? = 25.6, df= 1, p< 0.005), suggesting significant differences 
in genotype frequencies among examined populations. 
However, the allele-sharing genetic distance matrix failed to 
generate geographic clustering of individuals, and multi- 
dimensional scaling showed a poor fit of data into two 
dimensions (Stress = 0.30, R? = 0.53). 


DISCUSSION 


Genetic Differentiation of Argalis in Mongolia: Due 
to a lack of apparent topographic boundaries and yet highly 
variable morphology across populations, controversy 
surrounds the level of genetic differentiation among 
Mongolian argali populations (Allen 1940; Sopin 1982; Valdez 
1982; Geist 1991). This controversy has continued in large 
part due to the difficulties of obtaining genetic samples from 
the remote range of this species, and these difficulties also 
limited the sample size of our own analyses. Nevertheless, 
despite the limited sample size, our analyses of both mtDNA 
control region and nuclear microsatellites revealed significant 
genetic differentiation among sampled argali populations in 
Mongolia. 

The mtDNA control region phylogeny revealed two 
major groups — one composed of individuals from Altay and 
the other comprised primarily of individuals from Hangay/ 
Gobi — with an average sequence divergence (4.39%) similar 
to that observed between subspecies in other large mammals 
(Douzery and Randi 1997; Wooding and Ward 1997; Arctander 
et al. 1999; Matsuhashi er al. 1999), and greater than that 
typically seen among populations of Bighorn Sheep (Ovis 
canadensis; Ramey 1995; Luikart and Allendorf 1996; Boyce 
et al. 1999). This is also consistent with the speculation that 
habitat differences and a subtle geographic barrier (the 
Alakhnur Depression) have led to a reproductive isolation 
between Argali Sheep living in the Altay Mountains and those 
living in the Gobi desert (Sopin 1982). Interestingly, one Altay 
haplotype grouped with, but was basal to, the Gobi/Hangay 
group. Though this might be due to a low level of gene flow 
between the regions, incomplete lineage sorting seems a more 
plausible explanation. 


41 


GENETICS OF MONGOLIAN ARGALI SHEEP 


Within the Hangay/Gobi group, all Hangay haplotypes 
formed one monophyletic subgroup and the single Gobi 
haplotype was positioned outside this Hangay group. Although 
this suggests some differentiation between Hangay and Gobi, 
additional samples from Gobi are required to further address 
the phylogenetic relationship between these populations. 

Our nuclear microsatellite results also indicated 
significant genetic differentiation among the sampled 
populations. Although our sample sizes were small, the 
observed F , value of 0.056 1s similar to that found among 
natural populations of other large mammals (Roy ef al. 1994; 
Forbes and Hogg 1999; Paetkau et al. 1999; Gutiérrez-Espeleta 
et al. 2000). Moreover, although the individual allele-sharing 
genetic distance matrix did not generate meaningful geographic 
groupings (probably due to small sample size) the high 
percentages of “correct” assignments yielded in the likelihood 
assignment test suggest that allele frequency distributions 
differ among the three populations. Furthermore, the similar 
pair-wise F „and genetic distances and the even distribution of 
unassigned individuals suggest that these three populations 
are approximately equally differentiated from each other. 

Implications for Argali Taxonomy: Currently, two 
subspecies of argali are commonly recognized in Mongolia 
(Sopin 1982; Valdez 1982; Geist 1991): O.a. ammon (Altay 
mountain region) and O.a. darwini (Gobi desert region). 
These subspecific designations are supported by 
morphometric analyses (Kapitanova et al. 2004), and our 
mitochondrial control region phylogeny supports the 
distinction between argali from these regions. In contrast, 
although Hangay argalis are currently classified as 
O.a. ammon based on morphological similarities (Sopin 1982; 
Geist 1991), our mitochondrial analyses instead suggests that 
Hangay argalis are more closely related to O.a. darwini than 
to O.a. ammon (Tserenbataa et al. 2004). 

There are two possible explanations for this 
discrepancy. First, the more ammon-like morphology of 
Hangay argali may be due to the higher habitat productivity 
of the Hangay region relative to the arid Gobi desert. Second, 
because the mtDNA phylogeny only represents maternal 
descent, Hangay argali may be a hybrid form resulting from 
matings between large-bodied ammon males and small-bodied 
darwini females. This ‘hybrid origin’ hypothesis also can 
explain the approximately equal genetic distances that we 
obtained from nuclear microsatellite data. One way to test 
this hypothesis 1s to use a Y-linked marker to reconstruct the 
paternal lineage of Mongolian argalis. Moreover, since we 
were able to obtain only one sequence from Gobi argali, it is 
possible that the Hangay and Gobi populations represent two 
distinct subspecies of argali; this possibility requires further 

> with additional haplotypes. 


Implications for Conservation Management: Our 
mtDNA results showed that Mongolian argali haplotypes can 
be divided into two reciprocally monophyletic groups — one 
consisting of haplotypes found only in the Altay Mountains, and 
the other consisting almost exclusively of haplotypes from the 
Hangay Mountains and eastern Gobi desert. Due to the presence 
of one Altay haplotype in the Hangay/Gobi clade, argali in these 
two regions are not strictly reciprocally monophyletic, and 
therefore do not fit the definition of Evolutionary Significant 
Units (ESU’s) suggested by Moritz (1994). Nevertheless, our 
mtDNA results suggest two distinct, independent lineages, and 
the ESU criterion suggested by Moritz (1994) has been criticized 
for being overly stringent (Crandall et al. 2000; Fraser and 
Bernatchez 2001). Therefore, we tentatively recommend that 
argali in Altay and Hangay/Gobi be treated as two separate ESU’s 
for conservation purposes. 

Our microsatellite analyses indicated significant nuclear 
genetic differentiation among all three regions of Mongolia. 
Our mitochondrial control region analyses also showed 
differentiation between our single Gobi haplotype and all 
Hangay haplotypes, with the latter forming a single 
monophyletic clade. These results suggest that each area 
should be treated as a separate management unit (Moritz 1994) 
for conservation purposes. However, this recommendation 
should be considered tentative because sample sizes and areas 
surveyed were limited in this study (particularly for 
mitochondrial analyses), and because microsatellites can 
sometimes show significant differentiation across populations 
that may not be biologically meaningful (Hedrick 1999). 

Our results are mostly consistent with the mitochondrial 
NDS5 analyses of Tserenbataa et al. (2004), who found 
significant genetic differentiation between the Altay and 
Hangay/Gobi regions. However, Tserenbataa et al. (2004) 
found that the differentiation between these two groups was 
relatively weak, that there was no significant differentiation 
between Hangay and Gobi populations, and that haplotypes 
from any region did not form a monophyletic group. 
Tserenbataa eft al. (2004) concluded that there has been 
significant historical gene flow among the three regions of 
Mongolia and nearby areas of China and Russia, and that the 
entire region should be treated as a single ESU/subspecies. 
Our results indicate that the lack of resolution in the ND5 data 
of Tserenbataa et al. (2004) likely is due to the slow 
evolutionary rate of ND5 (compared to the control region and 
nuclear microsatellites) and incomplete lineage sorting rather 
than to the movement of individuals between regions. 
Nevertheless, studies that combine large sample sizes (as in 
Tserenbataa et al. 2004) with more sensitive genetic markers 
(as in this study) are needed before making firm conclusions 
about conservation units for Mongolian and other argali. 
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